﻿// OCLVectorMatrix.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include "OpenCLKernel.h"
#include <math.h>
#include <time.h>


const char* _kernel_vector_matrix_ = KERNEL(
	__kernel void vector_matrix(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
	__local float4 Asub[16];
	__local float4 Bsub[16][4];
	const int li0 = get_local_id(0);
	const int gr0 = get_group_id(0);
	const int nTiles = colA / (16 * 4);
	float4 acc = (float4)0.0f;
	int gk, k;

	for (gk = 0; gk < nTiles; gk++)
	{
		for (k = 0; k < 16; k++)
		{
			Asub[k] = vload4((gk * 16 * 4 + k * 4) / 4, inputA);
			Bsub[k][0] = vload4(((gk * 16 * 4 + k * 4) * colB + gr0 * 16 * 4 + li0 * 4) / 4, inputB);
			Bsub[k][1] = vload4(((gk * 16 * 4 + k * 4 + 1) * colB + gr0 * 16 * 4 + li0 * 4) / 4, inputB);
			Bsub[k][2] = vload4(((gk * 16 * 4 + k * 4 + 2) * colB + gr0 * 16 * 4 + li0 * 4) / 4, inputB);
			Bsub[k][3] = vload4(((gk * 16 * 4 + k * 4 + 3) * colB + gr0 * 16 * 4 + li0 * 4) / 4, inputB);
			acc += Bsub[k][0] * Asub[k].x + Bsub[k][1] * Asub[k].y + Bsub[k][2] * Asub[k].z + Bsub[k][3] * Asub[k].w;
			/*acc.x += dot(Asub[k], (float4)(Bsub[k][0].x, Bsub[k][1].x, Bsub[k][2].x, Bsub[k][3].x));
			acc.y += dot(Asub[k], (float4)(Bsub[k][0].y, Bsub[k][1].y, Bsub[k][2].y, Bsub[k][3].y));
			acc.z += dot(Asub[k], (float4)(Bsub[k][0].z, Bsub[k][1].z, Bsub[k][2].z, Bsub[k][3].z));
			acc.w += dot(Asub[k], (float4)(Bsub[k][0].w, Bsub[k][1].w, Bsub[k][2].w, Bsub[k][3].w));*/
		}
	}
	vstore4(acc, (gr0 * 16 * 4 + li0 * 4) / 4, outputC);
	return;
}
);


const char* _kernel_vector_matrix2_ = KERNEL(
	__kernel void vector_matrix2(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
	__local float4 Asub, Bsub[4];
	const int li0 = get_local_id(0);
	const int gr0 = get_group_id(0);
	const int nTiles = colA / 4;
	float4 acc = (float4)0.0f;
	int gk, k;

	for (gk = 0; gk < nTiles; gk++)
	{
		Asub = vload4(gk, inputA);
		for (k = 0; k < 4; k++)
			Bsub[k] = vload4(((gk * 4 + k) * colB + gr0 * 16 * 4 + li0 * 4) / 4, inputB);
		acc += Bsub[0] * Asub.x + Bsub[1] * Asub.y + Bsub[2] * Asub.z + Bsub[3] * Asub.w;
	}
	vstore4(acc, (gr0 * 16 * 4 + li0 * 4) / 4, outputC);
	return;
}
);


const char* _kernel_vector_matrix3_ = KERNEL(
	__kernel void vector_matrix3(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
	__local float4 Asub[16], Bsub[16];
	const int li0 = get_local_id(0);
	const int gr0 = get_group_id(0);
	const int nTiles = colA / (16 * 4);
	float acc = 0.0f;
	int gk, k;

	for (gk = 0; gk < nTiles; gk++)
	{
		for (k = 0; k < 16; k++)
		{
			Asub[k] = vload4((gk * 16 * 4 + k * 4) / 4, inputA);
			Bsub[k] = (float4)(inputB[(gk * 16 * 4 + k * 4) * colB + gr0 * 16 + li0],
				inputB[(gk * 16 * 4 + k * 4 + 1) * colB + gr0 * 16 + li0],
				inputB[(gk * 16 * 4 + k * 4 + 2) * colB + gr0 * 16 + li0],
				inputB[(gk * 16 * 4 + k * 4 + 3) * colB + gr0 * 16 + li0]);
			acc += dot(Asub[k], Bsub[k]);
		}
	}
	outputC[gr0 * 16 + li0] = acc;
	return;
}
);


// vector: colA; matrix: colA * colB
const int colA = 57600, colB = 2304;


bool floatEq(const float a, const float b)// 相等返回 1
{
	if (b == 0)
		return fabs(a) < 0.001;
	return fabs(1.0 * a / b - 1) < 0.001;
}


int main()
{
    float* A = (float*)malloc(sizeof(float) * colA);
    float* B = (float*)malloc(sizeof(float) * colA * colB);
    float* C = (float*)malloc(sizeof(float) * colB);
    srand(2);
    for (int i = 0; i < colA; A[i] = rand() & 0xF, i++);
    for (int i = 0; i < colA * colB; B[i] = rand() & 0xF, i++);
    for (int i = 0; i < colB; C[i] = 0, i++);

    cl_int status;

    OpenCLKernel* m_pOCL = new OpenCLKernel;
    if (0 != m_pOCL->OpenCLInitialize())
    {
        delete m_pOCL;
        m_pOCL = nullptr;
        return EXIT_FAILURE;
    }

    cl_command_queue queue = clCreateCommandQueue(m_pOCL->GetContext(), m_pOCL->GetDevices()[0], NULL, &status);
    GET_CL_ERROR(status, m_pOCL);

    //编译
    size_t nSourceSize = NULL;
    cl_program program = NULL;
    cl_kernel kernel = NULL;

	nSourceSize = { strlen(_kernel_vector_matrix_) };
	program = clCreateProgramWithSource(m_pOCL->GetContext(), 1, &_kernel_vector_matrix_, &nSourceSize, &status);
	GET_CL_ERROR(status, m_pOCL);
	status = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
	GET_CL_ERROR(status, m_pOCL);
	kernel = clCreateKernel(program, "vector_matrix", &status);
	GET_CL_ERROR(status, m_pOCL);

	//读内存
	cl_mem bufferA = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * colA, A, &status);
	GET_CL_ERROR(status, m_pOCL);
	cl_mem bufferB = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * colA * colB, B, &status);
	GET_CL_ERROR(status, m_pOCL);
	cl_mem bufferC = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_WRITE_ONLY, sizeof(float) * colB, NULL, &status);
	GET_CL_ERROR(status, m_pOCL);

	clock_t time_start = clock();

	//参数赋值
	status = clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufferA);
	GET_CL_ERROR(status, m_pOCL);
	status = clSetKernelArg(kernel, 1, sizeof(cl_mem), &bufferB);
	GET_CL_ERROR(status, m_pOCL);
	status = clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufferC);
	GET_CL_ERROR(status, m_pOCL);
	status = clSetKernelArg(kernel, 3, sizeof(int), &colA);
	GET_CL_ERROR(status, m_pOCL);
	status = clSetKernelArg(kernel, 4, sizeof(int), &colB);
	GET_CL_ERROR(status, m_pOCL);

	//工作组
	size_t* globalSize = NULL;
	size_t* localSize = NULL;
	cl_uint work_dim = 0;
	work_dim = 1;
	globalSize = new size_t[1];
	localSize = new size_t[1];
	globalSize[0] = colB / 4;
	localSize[0] = TILE_DIM;

	//执行
	status = clEnqueueNDRangeKernel(queue, kernel, work_dim, NULL, globalSize, localSize, 0, NULL, NULL);
	GET_CL_ERROR(status, m_pOCL);
	status = clFinish(queue);
	GET_CL_ERROR(status, m_pOCL);

	clock_t time_end = clock();
	printf("\nTime kernel : %d ms\n", time_end - time_start);

	//写内存
	status = clEnqueueReadBuffer(queue, bufferC, CL_TRUE, 0, sizeof(float) * colB, C, 0, NULL, NULL);
	GET_CL_ERROR(status, m_pOCL);

	// 返回并检查结果
	printf("Checking...\n");
	int correct = 1;
	for (int j = 0; j < colB; j++)
	{
		float tempSum = 0.0f;
		for (int i = 0; i < colA; i++)
			tempSum += A[i] * B[i * colB + j];
		if (!floatEq(tempSum, C[j]))
		{
			printf("Error at [%d], calculation: %f, reference: %f\n", j, C[j], tempSum);
			correct = 0;
		}
	}

	//结束并释放
	if (globalSize != NULL)
		delete[] globalSize;
	if (localSize != NULL)
		delete[] localSize;
	clReleaseMemObject(bufferA);
	clReleaseMemObject(bufferB);
	clReleaseMemObject(bufferC);
	clReleaseKernel(kernel);
	clReleaseProgram(program);
	clReleaseCommandQueue(queue);
	delete m_pOCL;
	free(A);
	free(B);
	free(C);
    return 0;
}

